Immigration as the main driver of population dynamics in a cryptic cetacean

Abstract Empirical evidence about the role and interaction of immigration with local demographic processes in shaping population dynamics is still scarce. This knowledge gap limits our capability to derive a conceptual framework that can be used to inform conservation actions. Populations exposed to nonstationary environment do not converge to a stable stage distribution, implying the need for evaluating the demographic role of both vital rates and stage distribution using appropriate tools. This is particularly important for species with larger generation times like cetaceans. We explored the relative demographic role of vital rates and population structure of a poorly known cetacean, the Mediterranean Cuvier's beaked whale, while accounting for the exposure to nonstationary environments. We performed a retrospective analysis through transient life table response experiments (tLTRE) using demographic rates and population structure of both sexes obtained from an integrated population model. The contribution of immigration to variation in realized population growth rates was 4.2, 7.6, and 12.7 times larger than that of female apparent survival, proportional abundance of breeding females with a 2‐year‐old calf, and proportional abundance of breeding females with a 3‐year‐old calf, respectively. Immigration rate and proportional abundance of breeding females with a 2‐ or 3‐year‐old calf explained, respectively, 65% and 20% of total temporal variability in realized population growth rates. Changes in realized population growth rate between successive years were mainly driven by changes in immigration and population structure, specifically the proportional abundance of breeding females with a 2‐year‐old calf. Our study provides insight into the demographic processes that affect population dynamics and in a cryptic cetacean. We presented an analytical approach for maximizing the use of available data through the integration of multiple sources of information for individuals of different distinctiveness levels.


| INTRODUC TI ON
Understanding the demographic consequences of changes in vital rates and population structure is important to assess factors that can threaten the viability of animal populations and to improve conservation actions (Morris & Doak, 2002). However, information on birth and survival is available only for 1.3% of the tetrapod species, and for 65% of threatened tetrapods we have no information about their vital rates (Conde et al., 2019). The sensitivity of population growth to changes in vital rates is affected by species' life history, with long-lived species that mature late and produce few or a single offspring, like cetaceans, that are expected to invest more in survival than in reproduction (Brault & Caswell, 1993;Saether & Bakke, 2000;Young & Keith, 2011). Despite many studies have shown a greater proportional sensitivity of population growth rates to adult survival in long-lived species, this has been questioned considering that other vital rates, like juvenile/immature survival and reproduction, can show larger temporal variability that may translate into a greater impact on population dynamics (Brault & Caswell, 1993;Manlik et al., 2016;Mills & Lindberg, 2002;Sergio et al., 2011).
Obtaining reliable estimates of demographic rates is critical for evaluating population dynamics and viability but, even when individual longitudinal data to estimate survival are available, the collection of additional data on fecundity rates and immigration can be difficult in hard-to-study species. The overall increased availability of data and analytical tools in the past decade has improved our ability to disentangle factors affecting the dynamics of animal populations.
Among these factors, immigration is notoriously difficult to record and quantify in wild animal populations (Williams et al., 2002).
Quantitative estimates of immigration rates have been made available only recently for some bird and mammal species, for which this vital rate has been confirmed as a key demographic parameter (Millon et al., 2019). While dispersal processes, emigration and immigration, have been recognized as important drivers of metapopulation persistence and local population dynamics (Hanski, 1999), empirical evidence about the role and interaction of immigration with local demographic processes in shaping population dynamics is still scarce, even for the most studied avian and mammalian populations. This empirical knowledge gap limits our capability to derive a conceptual framework that can be used to inform conservation and population ecology (Millon et al., 2019).
Immigration can be estimated through genetic analyses (Rannala & Mountain, 1997;Wang & Whitlock, 2003) or, more precisely, through demographic methods using data on marked animals at multiple populations or through a joint analysis of multiple data types.
The latter is performed within an integrated population modeling framework (IPM; Abadi, Gimenez, Ullrich, et al., 2010;Besbeas et al., 2002) that incorporates a common likelihood typically for population counts, fecundity, and encounter-reencounter data. In this case, information about immigration and its temporal variability contained in population count data can be derived because the other data sources contain information about the remaining vital rates, and because the IPM framework ensures an adequate representation of errors of the immigration parameter (Millon et al., 2019).

The joint use of estimates of demographic parameters from
IPMs and transient life table response experiments (tLTRE; Koons et al., 2016Koons et al., , 2017 has improved in the past few years our understanding about the relative demographic role of vital rates and population structure while accounting for the fact that populations are exposed to changing conditions, that is, nonstationary environments. Empirical evidence derived from the latter framework is, however, limited to a few species and case studies, mainly regarding bird populations (Koons et al., 2017;Nater et al., 2021;Schaub & Ullrich, 2021;Tenan et al., 2021). To our knowledge, no studies on marine mammals employed tLTREs to explore the demographic contribution of vital rates and population structure under nonstationary environmental conditions. A practical limitation of applying tLTREs is in fact the need for estimates of both vital rates and population structure, the latter being difficult to infer for many vertebrates. Here, we used multiple data types integrated into an IPM and tLTREs to study the population dynamics of an endangered cetacean, the Cuvier's beaked whale (Ziphius cavirostris). Beaked whales are the least understood among large mammals at a global scale and are particularly vulnerable to anthropogenic noise disturbances caused by geological and seismic surveying, military sonar, and naval traffic that can cause mass strandings (Li & Rosso, 2021). The paucity or lack of baseline data on distribution, population size and structure, and life history of beaked whales limits or impedes quantifying anthropogenic impacts at the population level and design evidence-based conservation plans (Hooker et al., 2019). The Cuvier's beaked whale is, however, one of the few beaked whale species for which long-term encounter-reencounter data from photo-identification are available and referred to some populations (Curtis et al., 2021;Hooker et al., 2019). The joint analysis of encounter-reencounter data and population count data allowed us to estimate latent parameters, like immigration, derive multiple fecundity metrics, and explore the relationship between vital rates and

K E Y W O R D S
Cuvier's beaked whale, fecundity, integrated population model, interbirth interval, nonstationary, Ziphius cavirostris

T A X O N O M Y C L A S S I F I C A T I O N
Demography population structure and the growth rate of a Mediterranean Cuvier's beaked whale population.

| Study species and data collection
The Cuvier's beaked whale is a medium-sized toothed whale with cosmopolitan distribution, and the only beaked whale species commonly found in the Mediterranean Sea. Cuvier's beaked whale has the deepest (2992 m; Schorr et al., 2014) and longest (222 min; Quick et al., 2020) recorded foraging dives among marine mammals. In the Mediterranean, Cuvier's beaked whales occur in deep waters (>200 m) and are often found over the continental slope (Moulins et al., 2007) where they feed primarily on mesopelagic cephalopods (Santos et al., 2001). Genetic analysis has indicated a high degree of differentiation between the Atlantic and the Mediterranean population (Dalebout et al., 2005). The area investigated along the years represents a core area for the species in the Mediterranean sea (Moulins et al., 2008;Tepsich et al., 2014) and it extended beyond the 1000-m isobath, covering an approximate area of 2000 km 2 . The surveys were conducted in Beaufort sea state ≤4 and swell <1 m. For each sighted group the geographical position was recorded and the group was described in situ on the basis of the minimum, the maximum, and the best estimation of the number of individuals. Photographic captures of whales were collected using auto-focus digital cameras equipped with a 70-300 mm stabilized zoom lens. The sampling protocol was to photograph randomly as many animals as possible in the group, from both right and left sides, whether or not photographs had already been taken of a particular individual. Whenever possible, photographic sequences were taken in order to capture the whole exposed flank of the whales, from the snout to the tail stock (Coomber et al., 2016).

Encounter-reencounter information of all photographed individuals
collected during the 16-year period derived from a total of 233 surveys (mean 14.6, median 18, range 2-22 surveys per year) annually performed from April to the end of October, with 85% of records collected from June to September.

| Photo identification and data selection
Each photograph was assigned a quality rating (Q-value) from 1 to 6 (similar to Gowans & Whitehead, 2001), based on the focus, angle of the animal relative to the sensor plane, the proportion of the frames filled by the animal, and the exposure (for details see Rosso et al., 2011). The Q-value was independent of the amount of markings on the individual. Only captures events described by Q ≥ 4 photographs were analyzed in this study, since Q ≥ 4 images were of sufficient quality to be used in photographic encounter-reencounter studies in this species (Rosso et al., 2011).
Photographs of each whale-from each (photographic) capture event-were compared with photos from previous sightings.
Recognition of individuals was achieved through the observation of natural markings (Rosso et al., 2011). If a photograph matched an individual that was already known in the catalog, the photo was assigned the whale's identification number. If not matched, the individual was given a new identification number and added to the catalog. All the photos were analyzed by the last author, who had 15 years of experience in photo-identification of Cuvier's beaked whales. Photographic collections from left and right sides were maintained separately, although 93% of individuals were photographed on both flanks.
Individual distinctiveness was assigned according to the number and size of marks (i.e., shape, notches, and nicks) taking values of 0

| Demographic data
Two types of demographic data were modeled: encounterreencounter data and population count data. The encounterreencounter dataset included a total of 124 unique individuals (81 males and 43 females; Table S5 in Appendix S1) with 136 and 329 resightings for females and males, respectively (Tables S6 and S7 in Appendix S1); only one individual with unknown sex was discarded from the analysis. Population count data were composed by annual sex-and/or stage-specific counts of the total number of individuals  Tables S3 and S4 in Appendix S1). Information on fecundity was embedded in the encounter-reencounter data, through information of young-of-theyear/calf-mother association included in the defined event codes (see below).
The IPM framework allows the estimation of latent demographic parameters, for which no or very few explicit data are available, such as immigration (Abadi, Gimenez, Ullrich, et al., 2010;Schaub & Fletcher, 2015). Inference was based on a joint likelihood derived by the multiplication of the single-dataset likelihoods described below.

| Modeling of population count data
A sex-specific stage-structured stochastic population model, that describes changes in population size over time as a function of demographic rates, was used in a state-space modeling approach (De Valpine & Hastings, 2002), where the observation process that generated count data is conditional on the state process. The stages for females were juvenile, prebreeding, breeding (with a young-of-theyear, 1-, 2-, or 3-year-old calf), nonbreeding, and (juvenile) immigrant (see below). The stages for males were juvenile, subadult, not toothed adult, toothed adult, and immigrant. Changes in the sex-and stage-specific population sizes between year t and t + 1 were modeled with Binomial and Poisson processes (Equations S1-S45 and S48-S53 in Appendix S1). The number of immigrants in year t was N Imm, FM,t ∼ Pois t where ɷ t is the expected number of immigrants (including both sexes) in year t (Schaub & Fletcher, 2015). Immigration (ɷ t ) was modeled with random year effects and independently of other vital rates over time (Equation S47 in Appendix S1). The expected number of immigrants was not sex-specific to reduce bias in the estimation of this latent parameter, assuming even sex ratio among immigrants (Equations S13 and S35 in Appendix S1). In addition, immigration was assumed to occur in the juvenile stage. Sexand stage-specific probabilities were defined for female apparent survival (ɸ F,t , equal for all stages; Section S.1.1 and S.3 in Appendix S1) from year t to t + 1, young-of-the-year and calf survival ( By,F,t , conditional on mother survival F,t ), and the probability of remaining with the mother for a calf from year t to t + 1 (F t , conditional on survival By,F,t and fixed to 1 for young-of-the-year and 1-year old calves; Figure S1 in Appendix S1). Note that in the multievent model for encounter-reencounter data (see below) we modeled the apparent survival probability, Bc, F,t, for 1-and 2-year-old calves that remain with the mother up to the third year of life. However, apparent survival probability for 2-year-old calves may be biased since true mortality cannot be distinguished from calf emancipation from its mother (Couet et al., 2019). In the state-space model we thus assumed that true survival for a 2-year old calf was equal to By,t , i.e., survival of young-of-the-year and 1-year-old calves that cannot yet emancipate, which gives (conditional on mother survival) Bc,F,t = By,F,t F t , and probability of remaining with the mother F t = Bc,F,t By,F,t (see Appendix S1 for further details). Time-invariant transition probabilities from juvenile to prebreeding female ( JuvPB,F ) and from prebreeding to breeding female ( PbBy,F , i.e., probability of first reproduction) were also defined, along with the time-variant breeding probability for nonbreeding females (γ t ). The model also included the apparent survival probability for juvenile and subadult males from year t to t + 1 ( JuvSub,M,t ), the apparent survival probability for adult males (both not toothed and toothed; Ad,M,t ) and the stage transition probability ( M , equal among stages and time invariant, following the best model structure selected for the encounterreencounter data, see below).
The observation model described the relationship between annual counts of juveniles, nonreproductive individuals (i.e., immatures that have not yet recruited as breeders, excluding juveniles, and nonbreeders that have reproduced at least once in the past), breeding females, and adult toothed males with true sex-and stage-specific population sizes using Binomial processes with detection probabilities obtained from the multievent model (Equations S54-S57 in Appendix S1). Detection probabilities were shared between the two submodels given that the observation process for count data was the same as for photo-identification encounter-reencounter data.
This, however, implies the assumption that sex-and stage-specific detection probability does not change between individuals with distinctiveness levels 1 and 2.

| Modeling of encounter-reencounter data
Individual encounter histories were modeled using sex-specific multievent capture-recapture models (Pradel, 2005), in order to accommodate uncertainty in state assignment. We assumed eight states and seven types of mutually exclusive events that could be observed for females, and five states and six event types for males ( Table 2).
We also assumed that all female breeding states were ascertained with certainty.
Probabilities of survival for females ( F ), young-of-the-year ( By,F ), and calves ( Bc,F ), as well as breeding probability (γ) where modeled with random year effects and a variance-covariance matrix for temporal correlation between the parameters (Equations S62-S66 in Appendix S1). Similarly to female probabilities, male apparent survival probabilities ( JuvSub,M and Ad,M ) were also modeled with random year effects and with a temporal correlation structure (Equations S67-S69 in Appendix S1). Temporal random effects were used for modeling encounter probabilities of females and males, in addition to the effect of year-specific sampling effort in both sexes and a stage-specific effect on male probabilities (Equations S70 and S71 in Appendix S1). For further details about the multievent model structure and selection see Appendix S1.

TA B L E 2
States and codes of mutually exclusive events that could be observed implemented in the multievent model for encounterreencounter data of the Cuvier's beaked whale population.

| Fecundity
Fecundity metrics embedded in the IPM are the time-variant breeding probability for a nonbreeding female (γ t ), the interbirth interval derived as 1 ∕̂ + 1, where ̂ is the average breeding probability, the probability of first reproduction ( PbBy, F, i.e., the time-invariant probability of moving from the prebreeding to the breeding female state), and the weaning probability (i.e., the complement, 1−F t , of the probability of remaining with the mother for a calf of 2-3 years of age, from year t to t + 1, conditional on calf survival By, F,t). In addition to these quantities, another two were derived taking advantage of the integrated population modeling framework. The reproductive rate (RR t ), expressed as the proportion of females within the study area that have a young-of-the-year among those that have already recruited as breeder (i.e., have already reproduced at least once): In addition, the weaning rate per female (WR t ) was derived from the breeding female apparent survival probability (ɸ F ), the young-of- the-year/calf survival ( By, F; conditional on breeding female apparent survival), and the probability for a calf of remaining with the mother (F) until the third year of life, considering the two possible scenarios for a calf to be weaned between 2 and 3 years of age (first row of Equation 2) or at 3 years of age (second row of Equation 2):

| Bayesian inference
The integrated population model was fit using a Bayesian formulation with Markov chain Monte Carlo simulations. Vague prior distributions were used for all parameters (see supplementary model code for details on prior specification in Appendix S1). Summaries of the posterior distribution were calculated from 80,000 posterior samples (burn-in = 5,000,000 iterations). We assessed convergence using the R diagnostics (Brooks & Gelman, 1998) that was <1.02 for all parameters. Posterior predictive checks were performed to assess the goodness-of-fit of the different (sub)models (see Appendix S1).

| Demographic influence on population growth rate
We performed a retrospective analysis to understand population dynamics in relation to changes in demographic rates of both sexes where ‖ ‖ denotes the sum of absolute values of vector elements (i.e., the 1-norm), A t is the Leslie matrix and N t is the population vector in year t. For details about the product of the Leslie matrix and the population vector see Appendix S1 (Section S.2).
The first tLTRE measures the contribution of temporal variability in demographic rates and population structure to the temporal variance of realized population growth rate λ t . If we define vector q t whose elements are the demographic rates and the proportional stage-structured population sizes (Ñ), where proportional means that they sum to one for each year, the contribution of variability in each demographic rate θ i,t and each component of structured abundance θ j,t to the temporal variance of λ t is: where sensitivity of λ t to change in each underlying demographic pa-  Table S1 in Appendix S2). Conditional on mother survival, mean young-of-the-year, and calf survival was 0.835 (0.602-0.973) with a temporal variability considerably larger than that of survival of weaned individuals (Figure 2a; Table S1 in Appendix S2). Temporal correlation between demographic rates was generally weak, and 95% CRIs for temporal correlation coefficients always encompassed zero, except for female apparent survival probability and young-of-the-year survival (Tables S2 and (1) (2) WR t = F,t−3 By,F,t−3 F,t−2 By,F,t−2 F,t−1 By,F,t−1 1 − F t−1 + F,t−3 By,F,t−3 F,t−2 By,F,t−2 F,t−1 By,F,t−1 F t−1 . (

| Demographic influence on temporal variance of population growth
Decomposing the variance of realized population growth rates using transient life table response experiments (tLTRE; Appendix S1, Equation 3) showed that immigration rate contributed most to temporal variability in realized population growth rates (65% of total variation), followed by female apparent survival (13%), proportional abundance of breeding females with a 2-year-old calf (12%), and proportional abundance of breeding females with a 3-year-old calf (8%).
All other parameters contributed less than 3% of total variation each ( Table 3).

| Demographic influence on changes in population growth between successive years
The application of tLTRE to changes in realized population growth rate between successive years (Appendix S1, Equation 4), with sequential changes that ranged from −0.23 to 0.11, showing that when population growth rates changed substantially (Δ > 0.1 or Δ < −0.1) the dominant driver was immigration rate in all cases.
However, changes in population growth rate greater than | 0.1 | occurred only three times (2004-2005, 2005-2006, 2014-2015) during the study period. Considering a slightly lower absolute threshold value of changes in population growth rate (i.e. Δ > 0.08 or Δ < −0.08), which occurred six times during the study period, confirmed immigration rate as the dominant driver (67% of the cases) followed by the proportional abundance of breeding females with a 2-year-old calf (33%). All other demographic parameters and components of population structure contributed the least (Figure 3).

| DISCUSS ION
We Not accounting for imperfect detection of breeding females with young-of-the-year/calf may lead to underestimating fecundity (i.e.,

F I G U R E 2
Estimates of year-specific demographic rates (a-e) and population growth λ t (f) of the Cuvier's beaked whale population. In (c), the red dashed line indicates the mean from the random effects model, and the red dotted lines indicate the 95% credible intervals of the mean. The vertical lines indicate the 95% credible intervals for the annual estimates. The horizontal black dotted line in (f) indicates population stability. the longest birth intervals to be observed (Barlow & Clapham, 1997).
Hoyle and Maunder (2004) derived a fecundity rate per mature individual at equilibrium calculated as the inverse of the number of breeders per recruit at carrying capacity, but did not account for imperfect detection. The approach used by Arso Civil et al. (2017) to estimate fecundity rate, based on the probability of a female giving birth conditional on a previous birth, was able to take account for imperfect detection, individual and temporal heterogeneity in re-sightings. However, neither this nor the most commonly used methods could be used to explore temporal variation in reproductive TA B L E 3 Estimated sensitivities and elasticities of realized population growth rate to changes in the underlying sex-specific vital rates, overall immigration rate, and proportional population structure (i.e., stage-specific proportions of abundance, evaluated at temporal mean values), and transient life stable response experiment (tLTRE) contributions to temporal variability in realized population growth rate, for the Cuvier's beaked whale population in the period 2004-2019. used applies also to the immigration parameter itself, which can be expressed and modeled as a rate relative to previous population size or directly as the number of immigrating individuals (Schaub & Fletcher, 2015) like in our case. Immigration is a latent parameter, and we assumed it occurred at the juvenile stage and with even sex ratio among immigrants in order to reduce estimation bias and not further complicate the model due to the limited data available to inform the life-history parameter estimates. To our knowledge, estimates F I G U R E 3 Sequential changes in realized population growth rate for the Cuvier's beaked whale population, and contributions of the changes in demographic rates and population structure. The annual difference in realized population growth rate equals the sum across all contributions. Colors highlight the parameters whose changes contributed the most to at least one change in population growth during the study period; all the other parameters with minor contributions are grouped and reported in black.

| Demographic influence on population growth
We performed a retrospective analysis through transient life  (Schaub & Kéry, 2021). Populations that are exposed to a nonstationary environment remain in a permanent transient state with no chance to converge to a stable stage distribution (Hastings, 2004).
Therefore, populations that are not in a stable stage distribution depends not only on demographic rates but also on the actual stage distribution (Schaub & Kéry, 2021). In addition, in species with longer phases of transient dynamics, for example, species with larger generation times like cetaceans, the impact of population structure becomes stronger and the results of transient and classical LTRE will differ more (Koons et al., 2016). Transient LTREs should thus be preferred in demographic studies on cetaceans.
In our study population, the contribution of immigration to variation in realized population growth rates was 4.2, 7.6, and 12.7 times larger than female apparent survival, proportional abundance of breeding females with a 2-year-old calf, and proportional abundance of breeding females with a 3-year-old calf, respectively. Immigration rate (average 0.043) affected the most temporal variability in realized population growth rates (65% of total variation). Proportional abundance of breeding females with a 2-or 3-year-old calf explained 20% of temporal variation in realized population growth rates, suggesting that local recruitment is another important driver of population dynamics of this cetacean. Note that weaning probability for 2-year-old calves was on average 0.615 across the whole study period, suggesting that the majority of calves weaned before reaching their third year of life with the mother. Changes in realized population growth rate between successive years were mainly driven by changes in immigration and population structure, specifically the proportional abundance of breeding females with a 2-year-old calf.
The role of immigration as a driver of the dynamics of cetacean populations has never been quantified or considered before. Brault and Caswell (1993) used classical (i.e., not transient) LTRE and found that variance in growth rate was due to variance in adult reproductive output. Here, we also highlight the importance of dispersal through immigration processes.

| CON CLUS IONS
The IPM framework here used combines information from different datasets and allows the simultaneous estimation of vital rates, including latent parameters like immigration, and population structure, while accounting for sampling and process covariance among demographic rates. Estimation of vital rates is associated with estimation of stage-specific and total abundance and realized population growth rates, as opposed to asymptotic growth rates that require assumptions like stationary environment and a stable stage distribution. In order to exploit all available information, we made use of data for individuals that do not have permanent marks and cannot be identified across multiple years, by incorporating them into yearspecific population counts. These data were formally integrated with encounter-reencounter data from photo-identification of individuals with permanent marks classically used for the encounterreencounter analysis. With regard to the latter, our multievent model incorporates reproductive parameters following Couet et al. (2019) but could be further complicated for other species to account for variable litter size and uncertainty on the timing at offspring independence (Cubaynes et al., 2021). The model could also be extended to include spatial information, that is, spatial explicit encounterreencounter data, which might account for dispersal.
Our study provides insight into the immigration and reproductive processes that drive population dynamics of a poorly known cetacean. Quantifying vital rates and population structure, along with anthropogenic impacts and their population consequences are of primary concern to inform knowledge-based conservation actions (Hooker et al., 2019;Li & Rosso, 2021). Deep-diving species appear in fact particularly susceptible to acute effects of mid-frequency active sonar exposure, and Cuvier's beaked whales remain the most frequently observed species in sonar-associated stranding events worldwide (Curtis et al., 2021;Hooker et al., 2019). This poorly known, cryptic cetacean has thus become the focus of a wide range of studies that are collecting demographic data. We presented an analytical approach for maximizing the use of available data through the integration of multiple sources of information for individuals of different distinctiveness levels. By jointly using IPM estimates and tLTRE (Koons et al., 2016(Koons et al., , 2017, we have also shown how this framework allows for more insights into the relative demographic role of different vital rates and population structure in a cetacean, while accounting for nonstationary environments. Overall, given the large number of population studies on cetacean species that routinely record encounter-reencounter information, the approach presented here can be applied to other ecological systems to better understand the demographic contribution of lifehistory traits and population structure. This will represent a step forward to improve the assessment of how different human activities jointly influence the persistence of cetacean populations in addition to the increasingly recognized impacts on animal behavior (Nabe-Nielsen et al., 2018).

ACK N OWLED G M ENTS
We thank the large number of students, CETASMUS interns, volunteers, research staff, and collaborators who have helped to collect photos of Cuvier's beaked whales over the years. Data collection would not have been possible without the logistical support of Francesco Rebagliati, Alessandra Somà, Nicola Aurier, Albert Sturlese, Frazer Coomber, Maurizio Wurtz and Gianluca Bozzo. We thank Giacomo Tavecchia for his help with defining the multievent models. We also thank the anonymous referees for the constructive comments which helped to improve this paper.

FU N D I N G I N FO R M ATI O N
None.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data used in this paper are available on the Figshare